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(57) ABSTRACT 

A system and method for analyzing a medical image to 
determine whether an abnormality is present, for example, in 
digital mammograms, includes the application of a wavelet 
expansion to a raw image to obtain subspace images of 
varying resolution. At least one subspace image is selected 
that has a resolution commensurate with a desired predeter- 
mined detection resolution range. A functional form of a 
probability distribution function is determined for each 
selected subspace image, and an optimal statistical normal 
image region test is determined for each selected subspace 
image. A threshold level for the probability distribution 
function is established from the optimal statistical normal 
image region test for each selected subspace image. A region 
size comprising at least one sector is defined, and an output 
image is created that includes a combination of all regions 
for each selected subspace image. Each region has a first 
value when the region intensity level is above the threshold 
and a second value when the region intensity level is below 
the threshold. This permits the localization of a potential 
abnormality within the image. 
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NORMAL AND ABNORMAL TISSUE 
IDENTIFICATION SYSTEM AND METHOD 
FOR MEDICAL IMAGES SUCH AS DIGITAL 
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BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention relates to systems and methods for 
analyzing medical images, and, more particularly, to sys- 
tems and methods for analyzing digital mammograms. 15 

2. Description of Related Art 

Many computer-aided diagnosis (CAD) schemes have 
been devised for mammographic image analysis [1-27]. A 
general review of digital radiography has been given by Doi 2Q 
et al. [1]. Many of these methods are based on multireso- 
lution techniques. 

Work related to the use of various multiresolution meth- 
ods for investigating mammograms includes Refs. 3, 11, 12, 

19, 23, and 26. Dengler et al. [11] use a difference of two 2 5 
Gaussians for the detection filter, and the final detection is 
based on a global threshold. Valatx et al. [12] generate a 
smooth approximation of the image with a (3 -spline expan- 
sion and apply a mixed distribution based local thresholding 
technique to both the raw and approximated image; the 30 
output image is formed by subtracting the two thresholded 
images. A calcification segmentation method is developed 
by Qian et al. [3] using two-channel and multichannel 
wavelet transforms [19], based on subband selection and a 
rescaling (thresholding) technique for feature detection [24]. 35 
Strickland and Hann [23] apply the wavelet transform at full 
resolution (no downsampling) and detect independently in 
two sets (HH and LH+HL) of three full resolution subband 
images. The detection results are combined, further 
processed, and the inverse wavelet transform is imple- 40 
mented. De Vore et al. [26] implement the standard wavelet 
transform, select the important subbands, and invert the 
transform after wavelet coefficient suppression. The result- 
ing image is empirically thresholded in order to remove the 
remaining background information. 45 

Various statistical approaches have been used to study 
mammograms [12-14, 18,21,23,27]. Wavelet domain coef- 
ficient probability modeling has also been utilized in other 
areas of research: selecting optimized coding methods [28, 
29], Gauss-Markov field representation [30-32], and texture 50 
identification [32]. 

It is known that film grain noise in mammograms is signal 
dependent [33, 34]. Typically, the accepted noise field for 
radiographs results from three independent components: ( 1 ) 
spatial fluctuations in the number of x-ray quanta absorbed 55 
in the screen; ( 2 ) spatial fluctuations in the screen absorption 
associated with random structural inhomogeneities in the 
phosphor coating; and (3) spatial fluctuations in film sensi- 
tivity due to the silver halide random distribution per unit 
area in the emulsion [35]. Many CAD methods have found 60 
it essential to carefully treat the image noise with a prepro- 
cessing step [3, 15, 22, 27, 36]. 

SUMMARY OF THE INVENTION 

It is therefore an object of the present invention to provide 65 
a system and method for identifying normal and abnormal 
tissue in medical images such as mammograms. 


It is an additional object to provide such a system and 
method that permit significant time savings in reading 
clinically normal mammograms. 

It is a further object to provide such a system and method 
for providing a second opinion strategy. 

It is another object to provide such a system and method 
having sufficient performance to detect a predetermined 
portion of the normal images with a low probability of false 
negatives. 

It is yet an additional object to provide such a system and 
method for detecting calcifications. 

These objects and others are attained by the present 
invention, a system and method for identifying normal tissue 
in medical images. Here the term normal is intended to 
define an image that does not contain a suspicious area, an 
image aberration, or small image medium defects. As radi- 
ologists spends an enormous amount of time investigating 
images lacking abnormalities, the invention can save a great 
deal of valuable time. This system and method may also be 
considered a second-opinion strategy, since the image can be 
declared “normal” by the detection system and method, and 
then reviewed by a radiologist, and thus the image has been 
analyzed twice. 

In a preferred embodiment, the invention addresses the 
detection of microcalcifications in mammograms, with a 
performance of detecting 40-50% of the normal images with 
a low probability of false negatives. 

The invention comprises the use of a multiresolution 
statistical model for normal tissue. This model is then used 
to make comparisons with local image regions. If a small 
region deviates significantly from the global model, it is 
flagged as potentially suspicious; if a region is in agreement, 
it can be discarded. The systematic identification of abnor- 
mal regions can be regarded as a detection algorithm that can 
be tested and evaluated using a standard database. If no 
suspicious regions are located, an image lacking any pathol- 
ogy can be identified by the detection process. 

A fundamental distinction exists between the techniques 
of the present invention and the prior art. Herein a multi- 
resolution approach is used as a simplifying device for 
statistical modeling, in order to show that a multiresolution 
statistical analysis has the potential for simplifying what has 
previously been considered an intractable statistical prob- 
lem. Specifically, the statistical interpretation of a raw image 
is very difficult, but is reasonably simple when applied 
separately to various resolutions of the image, after the 
decomposition into independent subspace images. 

Specifically, the method of the present invention com- 
prises the step of applying a wavelet expansion to a raw 
image. The raw image, which is typically in electronic form, 
comprises an array of sectors (e.g., pixels), wherein each 
sector has an intensity level. The wavelet expansion is for 
obtaining a plurality of subspace images of varying resolu- 
tion. 

The next step comprises selecting at least one subspace 
image that has a resolution commensurate with a desired 
predetermined detection resolution range. For example, if it 
is desired to examine for the presence of a neoplasm or 
calcification having dimensions in a particular size range, 
one or more subspace images are selected that encompass 
that size range. 

Next is determined a functional form of a probability 
distribution function (pdf) for each selected subspace image 
and an optimal statistical normal image region test for each 
selected subspace image. There is a test statistic associated 
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with the normal image region test that has some pdf. From 
this pdf a threshold level is established for the probability 
distribution function from the optimal statistical normal 
image region test for each selected subspace image. Prefer- 
ably this step is accomplished with the use of a “test statistic, 
” which will be described in the following. 

Finally, an output image is created, such as in electronic 
and/or visualizable form. A region is defined as comprising 
at least one sector, typically a plurality of sectors, against 
which the threshold level is compared. The output image has 
a first value (e.g., “1”) for each region when the region 
intensity level is above the threshold and a second value 
(e.g., “0”) when the region intensity level is below the 
threshold. This image then permits the localization of a 
potential abnormality within the image. 

Preferably, this method is taken a step further, although 
this is not intended as a limitation. The further step com- 
prises determining for the presence of a plurality of above - 
threshold regions within a predetermined larger area. Such a 
plurality of above -threshold regions can be indicative of a 
likelihood of abnormality. 

There are two distinctions between the present invention 
and previously reported statistical approaches, including the 
noise processing approach. First, the statistical analysis is 
applied to independent subspace images. The accepted noise 
components and the signal are lumped together; the aggre- 
gate is considered as a random field. Second, the focus is on 
the identification of statistical properties of normal tissue at 
multiple resolutions. This tissue provides the most abundant 
signal in the image, since in some medical images, such as 
mammograms, most of the image information, even when 
abnormalities are present, is associated with disease -free 
tissue. Regions that contain statistical aberrations that devi- 
ate from the normal model are considered as abnormal 
regions. This yields a multiresolution model that allows for 
features of interest, in the case of mammograms, 
calcifications, to be considered as outliers or perturbations to 
the global statistic of the applicable subspace image(s). This 
is an image domain approach that illustrates the power of 
taking a different view of the same information. An alternate 
method is to work entirely in the wavelet domain. 

The features that characterize the invention, both as to 
organization and method of operation, together with further 
objects and advantages thereof, will be better understood 
from the following description used in conjunction with the 
accompanying drawing. It is to be expressly understood that 
the drawing is for the purpose of illustration and description 
and is not intended as a definition of the limits of the 
invention. These and other objects attained, and advantages 
offered, by the present invention will become more fully 
apparent as the description that now follows is read in 
conjunction with the accompanying drawing. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a raw image, size 1024x2048, with 35-^wrn digital 
resolution, rescaled for viewing. The ROI containing a 
benign calcification is outlined. 

FIG. 2 a The 256-256 ROI at fill resolution, magnified for 
viewing, and FIG. 2d the raw image histogram. The histo- 
gram is from the fill resolution 1024x2048 raw image. 

FIG. 2b is the low-resolution f 5 ROI and FIG. 2e the 
associated histogram. Comparison of these two histograms 
reveals that the irregular multimodal structure is contained 
in the background, and not much discernible information is 
present in this image. 

FIG. 2c The total image resulting from adding the five 
subspace images, d 1 -d 5 . Most of the image detail structure 
is contained here. 


4 

FIGS. 3 a-c include the fine to coarse detail 
representations, d 1 -d 3 . The empirical histograms; of FIGS. 
3 f-j (solid lines) are to be compared with the estimated pdfs 
(diamonds). Note that the histograms represent the absolute 
5 value image distribution. Points have been skipped in the 
theoretical plots to avoid overlap and confusion. 

FIG. 3 d-e are the fine to coarse detail representation. d 4 
and d 5 . The empirical histograms are to be compared with 
the estimated pdfs (diamonds). Note that the histograms 
10 represent the absolute value image distribution. Points have 
been skipped in the theoretical plots to avoid overlap and 
confusion. The benign calcification, noticeable in the d 4 
subimage, may be considered as an outlier to the normal 
tissue model at this resolution. 

5 FIGS. 4a and 4 b are the summary statistics for the d 3 and 
the d 4 subspace images, respectively. These cases are a fair 
representation of the theoretical and empirical agreement for 
all images studied. 

FIG. 5 is the detection flow chart for mammography. 

20 FIG. 6 is the raw image 2048x2048 pixels scaled by a 
factor of 2/5 for viewing purposes. The arrow points to the 
region containing the biopsy-proven cluster. 

FIG. 7 is the total combined detection mask projected into 
the sum of the first five detail images. 

25 FIG. 8 is the d 3 detected image, binary output 

FIG. 9 is the d 4 detected image, binary output. 

FIG. 10 shows a family of probability functions for 
N=l-4. Note that the curves spread out and become more 
bell shaped for larger N. 

30 FIG. 11 is an example of a case where N=1 for the first 
three detail images, top to bottom, respectively. The data 
represent about 5xl0 6 pixels, where the empirical histo- 
grams (solid) are to be compared with the theoretical esti- 
mate (dashed). The absolute value data are displayed. 

35 FIG. 12 is an example of a case where N=2 for the first 
three detail images, top to bottom, respectively. The data 
represent about 5xl0 6 pixels, where the empirical histo- 
grams (solid) are to be compared with the theoretical esti- 
mate (dashed). The absolute value data are displayed. 

FIG. 13 is an example of a case where N=3 for the first 
three detail images, top to bottom, respectively. The data 
represent about 5xl0 6 pixels, where the empirical histo- 
grams (solid) are to be compared with the theoretical esti- 
45 mate (dashed). The absolute value data are displayed. 

FIG. 14 is a generic test statistic curve indicating an 
arbitrary threshold. If moved to the right, the first type of 
error is reduced, and if moved to the left, the second error is 
reduced at the expense of the first. 

50 FIG. 15 illustrates generic test output for an image having 
no calcifications, where the squares enclose suspicious 
areas. This is an example of moving the threshold too far to 
the left. 

FIG. 16 illustrates detection output image with no calci- 
55 fications present, where the squares enclose suspicious 
areas. This is an example of moving the threshold to the 
right. Note that a relatively small amount of areas are 
returned. 

FIG. 17 shows a breast region excised from the off-image 
60 noise field. This makes comparisons between the multireso- 
lution properties of the breast field with the off-breast noise 
field. The detection method enables the breast region to be 
labeled; the white line marks the separation. The image has 
been overcontrasted to illustrate the concept. 

65 FIG. 18 is a large image section, with the arrow indicating 
a malignant calcification cluster and associated histogram, 
upper right. 
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FIG. 19 is the associated d 3 image and histogram (solid) 
and theoretical estimate (dashed), upper left. The absolute 
value data are used to create the histogram. Note the 
regularity of the histograms. 

FIG. 20 is the associated d 4 image and histogram (solid) 
and theoretical estimate (dashed), upper left. The absolute 
value data are used to create the histogram. 

FIG. 21a is the d 3 and FIG. 21 b d 4 test statistic curves: 
theoretical (dots) and empirical (solid). For this case a 
functional form was assumed. 

FIG. 22 shows the detection flow for one stage. 

FIG. 23 shows the detection output resulting from pro- 
jecting the binary mask into the sum of the first five detail 
images. This image is ready for further processing if desired. 

DETAILED DESCRIPTION OF THE 
PREFERRED EMBODIMENTS 

A description of the preferred embodiments of the present 
invention will now be presented with reference to FIGS. 
1-23. 

I. Image Information 

The images under investigation are film mammograms, in 
this particular embodiment digitized at 35-^m-per-pixel 
resolution with 12-b precision, using a DuPont NDT Scan II 
Film Digitizer, although this is not intended as a limitation. 
The use of large image sections avoids having the back- 
ground signal influence the model; this is preferable for 
reliable statistical analysis. The term “image” refers to these 
large sections. Any section is approximately the largest 
rectangle with power-of-two dimensions that can be 
inscribed within the breast boundary. These are manually 
excised in a particular embodiment. The larger images 
correspond roughly to sections 7x7 cm, and the smaller 
images correspond to 7x3.5 cm. Whole mammograms have 
also been tested. 

In an alternate embodiment (see FIG. 17), the tissue 
region can be electronically excised by assuming that two 
noise fields are present in the image: the tissue field that 
contains all the signal information within the tissue region 
and the off-tissue noise field that consists of two 
components, a generally random noise field and bright 
artifacts or anomalous regions. First the tissue region is 
separated from the random noise field, and then the largest 
remaining contiguous region containing signal is assumed to 
be the tissue region. Everything but this region is set to zero. 
A complete discussion of this technique is found in Ref. 51. 

The choice to work with high-resolution and high- 
dynamic-range images is prudent for two main reasons: (1) 
35-^m digital resolution closely matches the Nyquist sam- 
pling rate for the file used, and (2) a higher bit rate image 
more closely resembles the intensity distribution of the 
analog image. (If we err in the needed bit rate, it is better to 
err on the side of too much rather than too little.) 

A database utilized in the present investigation contains 
over 100 mammograms with combinations of normals (no 
pathologies) and abnormals (images with biopsy-proven 
calcification clusters), all with varying parenchymal 
densities, as described by Kallergi et al. [37]. Images from 
this database were used by Zheng et al. [25] at lower 
resolution (images reduced to 105-jum resolution) for calci- 
fication detection. We selected 28 at random and 2 specifi- 
cally. The 2 images not selected at random were deliberately 
picked because they contain very subtle clusters, and it is 
important to evaluate the detection performance with limit- 
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ing cases. Of the 30 images studied, 17 are clinically 
abnormal and 13 are pathology free. This means that 2 of 17 
(«12%) of the abnormal images in this study are difficult 
detection cases; this is well above the anticipated number of 
5 such cases likely to arise from a large database. 

II. Wavelet Expansion and Primary Statistic 
A. Multiresolution Expansion 

The multiresolution approach allows for the image expan- 
io sion to take the form of a sum of independent subspace 
images given by the identity [38] 

fo=(fo-fl)+(fl-f 2 > ■ ■ ■ ( 1 ) 

where [39, ch. 5] the f 2 image is the next coarser represen- 
tation of f 0 and f 2 is the next coarser representation of f 1 . 
Specifically, f 1 is a half resolution (lower half of the original 
frequency spectrum) version of f 0 and f 2 is a fourth- 
resolution (lower fourth of the original frequency spectrum) 
version of f 0 . The image that contains the difference in 
“information” between the successive images f 0 and f 1? is 
designated by d-^fo-f-L. This pattern is continued to yield 

f 0 =d 1 +d 2 + . . . +dj+fj. ( 2 ) 

25 The interpretation of this expression in terms of images is 
that the [fine detail-^coarse detail] of the image is contained 
in terms [d-^dy] and f '■ is a smoothed (blurred) version of f 0 . 
The important observation is that orthogonality of the sub- 
spaces ensures independence of the various images in the 
30 decomposition. This type of expansion and the relative 
frequency spectrum for each image are discussed by 
Daubechies [39, p. 332]. Each d,- image results from three 
independent components or subband images in the wavelet 
domain related to features oriented in the horizontal, 
35 vertical, and diagonal directions [39, pp. 313-320]. The 
image domain approach allows for the three directional 
subband components of a particular image resolution to be 
combined and observed simultaneously. 

The forward wavelet transform is performed with a sepa- 
40 rable kernel two-dimensional pyramid downsampling 
scheme. Each d- image is constructed from inverting the 
three appropriate subband wavelet domain images 
(horizontal, vertical, and diagonal components) of a particu- 
lar decomposition level j, with upsampling [39]. The for- 
45 ward and reverse transforms are executed with circular 
(periodic) boundaries. This wavelet decomposition and 
reconstruction method ensures that the dy images are linearly 
independent (lack redundancies) and ensures the image 
identity sum (2) be exact. The resulting dy image contains a 
50 limited bandpass (rectangular) section of the original fre- 
quency spectrum. For this analysis we use a symmlet 
wavelet with 12 coefficients. This wavelet is nearly sym- 
metric and is chosen because it roughly resembles the 
intensity profile of calcifications. This is supported by 
55 empirical observations and is only a general approximation. 
This does not imply that the symmlet functions are truly 
matched filters for every calcification profile, but merely 
indicates that this basis is probably a better choice than the 
Haar, D4, or other similar basis. 

60 For illustration purposes, we consider an image expansion 
with j=5. FIG. 1 presents a scaled version of a 1024x2048 
section of a mammogram with an arbitrary 256x256 region 
of interest (ROI) outlined. This region will be used for 
demonstration purposes, but the numerical analysis is per- 
65 formed on the entire fill section image. FIGS. 2 and 3 
contain the raw ROI, the low-resolution f 5 component, the 
sum of the five detail images, and the fine to coarse detail 
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information. Important observations stem from this decom- 
position include: (1) Image f 5 is a very blurred version of the 
raw image, and the image structure is hardly discernible (see 
FIGS. 2 a and 2b), (2) the total detail image (the sum of all 
the detail images) contains the important detail structure of 5 
the image (compare FIG. 2 a with 2c), (3) the information 
needed for calcification detection is contained in some of 
these images, mainly the d 4 image in this case (see FIG. 3d), 
and (4) the simplicity of the probability distributions for the 
detail images is apparent from the histograms of FIG. 3. 

B. The Family of Probability Distribution Functions 
1. Mammographic Applications 

Applicants have found that the expansion probability 
functions can be approximated by choosing the appropriate 
integer parameter N in the function: 

f ^_k exp(-JkW) ^ IN -l -2l(2\x\k) 1 
p (x) = 22N _ { 1 )! 

/=o 

where x represents pixel values, and k is a scaling parameter 
that must be determined, but is of little importance for 
practical purposes (see FIG. 10 for examples of N=l-4). The 
above equation results from application of an inverse Fou- 
rier transform of the function 


8 

(termed a difference image) for locating other than random 
changes (such as longitudinal studies in lung imaging for 
detecting nodules, in mammography for detecting tumors 
and architectural distortions, or manmade changes in natural 
scenes). This observation follows directly from Eq. (3) [47]. 
It would also seem plausible to relax the restriction on N in 
Eq. (3) and let it take on any positive value resulting in the 
most general probability model. 

C. Probability Model, N=1 

Based on our study of the 30 mammograms (more than 
5x30=150 d 7 - images) from the database, a good approxima- 
tion for the histogram associated with the individual du- 
nnage probability density functions (pdf) is given by 



where y ■ represents an arbitrary pixel value and c ■■■ is a 
constant. The index j is to emphasize that the constant c . is 
20 different for each d } - image. This distribution is known as a 
Laplace distribution [40] and is the primary statistic. If the 
random variable (rv) Y is distributed as hfy), then a new rv 
X=|Y| is distributed as a negative exponential distribution 
given by 
25 

p(xj; cj ) = — exp (—1; xj > 0 
c j V c j / 


In an initial mammography study [47, the disclosure of 
which is incorporated herein by reference] it was assumed 
that N=l. However, subsequent studies [48] indicate that this 
work was a special case of a more general formalism. The 
framework for this generalized expression is also an exten- 
sion of prior theoretical work [49]. The statistical analysis 
can be considered as the study of incremental change. 

The parameter N can be found by using the data from the 
expansion images by forming the empirical histograms. 
These data (empirical histogram) can be input into standard 
regression analysis that will find the best estimate of N in 
relation to it. [Regression analysis will find the best N that 
minimizes the error between the actual histogram curve and 
the theoretical estimate given in Eq. (3).] An observation 
stemming from previous investigations indicates that N=1 
most often, preceded by N=2, and then by N=3. Images and 
probability functions for N=l-3 are provided in FIGS. 
11-13, respectively. This inference is qualified by noting 
that it is based on the analysis of a limited data set. 

The statistical test in the most general case may change 
with N (but not necessarily). This implies that the likelihood 
function (see [47] for this development) must be generated 
on a per-image basis. One way to do this is to calculate the 
test for the anticipated range of N (perhaps N=l-5) and 
catalog the test procedure. Once N is determined, the test is 
known. The test used is based on the Neyman-Pearson 
lemma and is based on forming the ratio of two likelihood 
functions. That is, it is a test designed to disprove a null 
hypothesis that results in making a calculation (test statistic) 
in a small region of interest in the image. 

2. Extensions to Other Images 

The wavelet expansion images of complicated, predomi- 
nantly low-frequency, fields such as mammograms, chest 
radiographs, or natural scenes represent small changes as the 
image is blurred in a coarser version of itself (see Eqs. 1 and 
2 in [47]). This statistical analysis can be considered a study 
of incremental change. It would seem plausible that a similar 
analysis would translate to the study of the difference 
between images of the same scene acquired at different times 


This function is used for computational purposes and 
amounts to using the absolute value of the d y - image. This 
tacitly implies that the dj image has a symmetric pdf, and 
that the transform merely folds the left-of-zero part onto the 
right-of-zero part. The justification is that the d- image 
contains no intensity bias (no dc component); in general, 
there is no reason to expect positive pixel values to be 
favored over negative values. In Ref. [41], the same pdf was 
used as a distance measure for determining the local simi- 
larity between two similarly structured images for registra- 
tion purposes. The connection between the work in [41] and 
this analysis is suggested by the image identity (2). 

Following maximum likelihood arguments [42] for inde- 
pendent samples of the variable x J7 the parameter c y from (4) 
can be estimated by the average value of x- 

c j =<x j> 

We have assumed that the samples of x y - are independent; 
this is certainly not the case, but it does not pose a serious 
problem, as can be seen from Table 1. These results are for 
the image in FIG. 3 and are typical of results of all of the 
cases studied. The calculated averages of x - obtained directly 
from the data are in good agreement with the values obtained 
from the least-squares analysis of c ■. The empirical data are 
compared with the theoretical pdfs in FIGS. 3a-3e. 


TABLE 1 


PROBABILITY DISTRIBUTION PARAMETER c,- AND IMAGE 
AVERAGE (x ; ) FOR IMAGE dL- 

< 3 , 

di 

d 2 

d 4 

d 5 

Cj 

15.9 

26.6 32.6 

40.0 

59.8 

(*,-) 

15.8 

26.0 31.1 

38.1 

57.2 


We find empirically that the d 3 and d 4 images are most 
pertinent for calcification detection at this digital resolution. 
The wavelet transform and subimage selection are similar in 
some respects to other tried approaches [3], [19], [23], [26] 
in that the subband images (in the wavelet domain) are 
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selected a priori. For example, the transformation applica- 
tion is not like [23], but the selection is similar. Our method 
is based on using two independent images after wavelet 
inversion rather than combining the d 3 and d 4 components. 

By considering the size of the image compared to the 5 
number of pixels contained in a cluster, it follows that the 
cluster has a minimal effect on the global statistic. (There are 
roughly 5xl0 6 pixels in the image and about 2000 pixels for 
an average calcification cluster.) Therefore, the primary 
statistic can be considered as the model for normal tissue. 10 
This statistic is useful for developing robust statistical tests. 

D. Generalized Method for Determining the Operating 
Parameters 

As stated above, determining the functional form of the 
probability function for each expansion image determines 15 
the normal image region test and results in the test statistic, 
another probability curve that must be estimated. The test 
statistic is a calculation that is made in small image regions 
(for example, 8x8 or 16x16 pixels) across the applicable 
expansion images. In order to determine the detection 20 
threshold the test statistic distribution must be sampled in 
each expansion image. This implies making the calculation 
in various image regions (establishing a sample grid) and 
estimating the distribution. 

This test statistic distribution can be estimated in two 25 
ways: (1) assume a parametric form as in Section V, or (2) 
assume a nonparametric form and use a kernel density 
estimation approach [50]. The second method is more pow- 
erful if the data are not exact or if they differ from the 
assumed parametric function in the first case. The major 30 
point here is to estimate the normal tissue behavior of the 
test statistic; the technique used to estimate this distribution 
must “wash out” the presence of the calcifications in the 
curve (if present in the particular case). Both techniques will 
work; if appropriate, the former is more powerful if the data 35 
agree with the functional form, and the latter always will 
perform providing the bandwidth parameter in the kernel is 
selected properly. 

After estimating the test statistic curve, the operating 
thresholds must be determined. This is accomplished by 40 
making an initial guess and running an experiment on 
ground-truth data. These are data that are understood. The 
images are known either to be truly normal or to be 
abnormal, and if pathologies exist, the locations are known. 
This training data set must represent a fair cross section of 45 
images from the imaging environment in which the detec- 
tion method will eventually operate. The detection algorithm 
is implemented, and the results are counted. 

In detection problems of this kind there are two opposing 
errors: (1) Assume the area is abnormal when it is normal; 50 
or (2) assume the area is normal when it is not. The 
connection between the threshold and these error rates is 
easily assessed from the generic test statistic curve (FIG. 
14). This curve is an approximation of the normal test 
statistic behavior. In theory there is an additional curve that 55 
describes the abnormal behavior; we do not know this curve, 
but hypothetically it will overlap the test statistic curve in the 
right tail section, where more overlap implies more confu- 
sion in the decision process (in effect this procedure is 
assessing this unknown curve). The effect of moving the 60 
threshold can be observed in the sample detection output 
(FIGS. 15 and 16). As the threshold is moved to the left, 
more regions will be declared as abnormal that are really 
normal (this is in a direction where the hypothetical curve of 
the abnormal tissue is dying out; FIG. 15). As the threshold 65 
is moved to the right, all areas that are normal will be 
declared normal, and the infringement on the hypothetical 
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curve will result in declaring abnormal image regions as 
normal (FIG. 16). This approach allows the calculation of 
the first error; the portion of the area under the test statistic 
curve from the threshold to infinity (to the right) results in 
the fractional error of declaring an image area abnormal 
when it is normal. 

In order to find the optimal tradeoff between the two 
errors, the experiment is preferably run many times, and the 
detection results are counted. Plotting the results will indi- 
cate the best operating parameters. There are two ways that 
this can be accomplished: (1) vary the two thresholds (one 
for each expansion image) in unison; or (2) vary the thresh- 
olds independently. In addition, if the tests are different for 
different N, further subclassification may be necessary, 
where all like N are segregated into one class, and thresholds 
adjusted by the value of N independently (this is the most 
general case). 

Thus the method and system of the present invention 
accomplish the following: 

The technique uses a multiresolution analysis that allows 
that patient to become its own reference for normality 
for a given digital system and imaging protocol, per- 
mitting great generality. 

The multiresolution expansion supports the parametric 
analysis, which would otherwise not be possible. 

The actual image data dictate the test, as opposed to 
forcing a test on the problem. 

Reasonable estimates of the error rate can be made. 
Therefore, an assessment can be made in advance of the 
possibility of declaring an image region abnormal when 
it is normal. 

The approach comprises a statistical analysis of incre- 
mental change, which is a vital component of the 
analysis. 

The analysis of incremental change should translate into 
the serial analysis of difference images (the difference 
between the same scene acquired at different times), 
with applications in longitudinal (time series) 
mammography, chest radiography, magnetic 
resonance, and general natural scene analysis. 

E. Background Analysis 

Traditionally, the background has been loosely defined as 
the information that is responsible for the irregular contrast 
qualities and variegated appearance of the image. It has been 
recognized that the background information presents a prob- 
lem for detection techniques [12]; in particular, dense tissue 
(which appears bright) is a major factor in making certain 
mammograms harder for both machine and radiologists to 
interpret [43], [44]. When the low-frequency constituent is 
separated (this comes from discarding the f constituent; see 
FIG. 2b), the dense and fatty areas (fatty regions appear as 
dark areas in the raw image) appear with equal contrast; this 
is illustrated in FIG. 2c. A conclusion may be that the 
radiologist should be given the option to view the entire 
image in conjunction with any combination of detail images 
with or without the background subimage included. In 
essence, this method for removing the background is a linear 
technique that restores contrast to the entire image. 

Although our studies on modeling the low-frequency f 5 
and f 6 distributions have not yet yielded any firm analytical 
conclusions, some qualitative observations can be made as 
to general trends associated with many images. The irregular 
structure of the raw image pdf appears multimodal, as shown 
in FIG. 2a. Note that the irregular structure of the raw image 
distribution is also present in the low-frequency image pdf 
and can be observed by comparing FIGS. 2a and 2b. A 
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in Table 2 are expected due to the integer rounding; the 
distribution images are saved as short integers. 

An error of the first kind or FP rate can be estimated from 
this pdf prior to detection processing. Again, it is assumed 
that the calcified regions have a minimal effect on this 
distribution and can be considered as outliers located in the 
far right tail region. The FP rate (the fractional number of 
false calcifications per image) can be obtained by 

Pf= f f(c;a,j3)dc 


where x denotes the threshold. This equation deserves spe- 
cial consideration. The test criterion given by <x»y is 
obtained as follows: (1) select T=y, (2) pick a value for Py, 
and (3) solve this equation for t. The total expected number 
of FPs in the entire dy image can be approximated by 

.FP(total)-PyX (number of pixels in reduced image) 

It should be emphasized that this is an estimation that may 
be obtained as an average after processing many images. 
The intriguing aspect is that the FP rate Py can be set prior 
to detection. However, this does not completely specify the 
error prediction, since no preprocessing estimation can be 
made concerning the error of the second kind or FN rate. The 
FN rate follows from the calcification distribution, which is 
unknown. 

VI. Localized Normal Tissue Recognition 

The detection technique is implemented by shifting an 
8x8- or 16x1 6-pixel search window though the d 3 and d 4 
images, respectively. A detection flow diagram that illus- 
trates the various stages of processing is given in FIG. 5. The 
intent is to match the search window size to the average 
spatial extent of the calcifications that may exist in each 
subimage. When the spatial extent of the wavelet function 
and calcifications are similar, the response (in the d ; - image) 
is maximized and the area is flagged as suspicious. 

For early cancer detection calcifications with spatial 
extent less than ~0.5 mm are most important for clinical 
diagnosis. This corresponds to calcifications ranging 
roughly from 16 pixels to three or four pixels in diameter 
(«0.1 mm), and the search window is matched to this scale. 
We assume that calcifications smaller than this are not 
discernible. The window is shifted with a 50% overlap in 
both spatial dimensions during the search; this is to reduce 
the risk of missing a feature (calcification). When a region 
is assumed normal (accept the null hypothesis), it is set to 
zero. If the null hypothesis is rejected (accepted as 
suspicious), the region is left intact. This is how potentially 
small calcified regions are detected by default. Following 
the independent detection, the images are combined and the 
total detected image results as the output. In effect, the dual 
output combination can be viewed as a mask. From this 
mask, any d ; - image combination or even the raw image can 
be returned as the output image. This can be accomplished 
by making the total output image into a binary image (ones 
or zeros) and simply multiplying by the desired type of 
output. This is important if further processing is desired 
because the calcified regions can be returned with full 
resolution and detail. For demonstration purposes, we 
choose an image with a wide variety of aberrations illus- 
trated in FIG. 6. The arrow indicates a biopsy -proven 
calcification cluster. The detection output is from the sum of 
the first five d } - images. These results are illustrated in FIG. 
7. 
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The detection scheme takes into account that calcifica- 
tions have a spatial extent or connectivity quality; a calci- 
fication appears as a clump of large pixel values and normal 
regions have a diffuse distribution of pixel values. The 
5 regions considered are matched to spatial scale by adjusting 
the search window size accordingly. The window size and 
shift increment are a compromise. If the window is greater 
than half the size of the feature and the shift greater than half 
the window size there is a possibility that the feature will be 
10 missed. Assume that the smallest object to be detected has 
a spatial extent M, a 2Mx2M window must be used and 
shifted with an increment of M to ensure that the object will 
not be missed. The window is most sensitive when the 
feature fills it entirely, which is not generally expected here. 
15 An alternative method would be to scan the image with the 
limiting window size and single pixel increments, but this 
may add many FPs to the outcome. 

Although not apparent, there is redundancy built into this 
detection scheme. This can be assessed by looking at the d 3 
20 and d 4 detected images prior to recombination, see (FIGS. 8 
and 9). There are flagged regions in the d 3 image that are not 
flagged in the d 4 image and vice versa. This indicates that the 
wavelet response to the feature was stronger in the respec- 
tive image. However, in some regions the test is triggered in 
25 both images at roughly the same spatial location. This 
indicates that the wavelet function response is similar to 
both images (relative to the background and window size). 
Thus some calcified regions have the possibility of being 
detected in both images, which represents a redundancy. 
30 This can be viewed as a safety measure. 

The focus of this detection scheme is very localized. 
However, a possible sign of early cancer is the presence of 
a microcalcification cluster, and this is of more clinical 
concern than isolated events. A single cluster is defined 
35 roughly as 3-5 microcalcifications assembled within a 
square centimeter. This definition implies that one or two 
calcifications within a square cm are not clinically impor- 
tant. In order to reach the goal of recognizing images that are 
normal from the clinical point of view clearly requires 
40 another stage of processing. 

The additional stage of processing is needed to eliminate 
FP normal diffuse regions on the order of a square cm. The 
FP normal diffuse region may result from isolated events 
within a square cm proximity: two PF calcifications (flagged 
45 regions that are normal) and one true calcification (a cor- 
rectly flagged region); vice versa; or three FP regions. 

VII. Experiment, Analysis, and Evaluation 

In this section the experimental implementation of the 
50 ideas expressed in Sections II. D and V are shown. In 
addition, the error rates are demonstrated to be assessable in 
advance by knowledge of the test statistic. 

A. Experiment 

The localized normal region detection must be conducted 
55 such that the FN and FP rates are optimized. These rates are 
in opposition in that decreasing one causes an increase in the 
other. In terms of the threshold, if t is set low enough the FN 
rate can be reduced to zero, but then the FP rate is high. So 
the problem is to adjust the threshold. We want the threshold 
60 as high as possible while keeping the FN rate essentially 
zero. This optimum value can be found by probing the 
detection operating characteristics. We do this by processing 
the images 5 times, each time with a slightly higher 
threshold, or equivalently, a lower value for Py. 

65 B. Analysis Method 

The evaluation of the local area detection method for each 
of the trials was performed by a resident radiologist using 
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three figures of merit. [First, for clarification, a true positive 
(TP) with respect to an isolated calcification is defined as: a 
calcified region, benign or malignant, that has not been set 
to zero. The TP cluster follows from this definition also.] 
The figures of merit are: (1) the isolated FP calcifications per 
image; (2) the number of TP clusters; and (3) the number of 
FN clusters. The cluster analysis is based on the biopsy 
verified ground truth files, and the results are presented as 
averages. There are many methods used for counting clus- 
ters; consequently, the technique used here requires a brief 
explanation. 

Following from the definition of a cluster (as defined 
previously), if three events are located (this includes FPs or 
TPs) within a square cm, the region is classified as a cluster. 
If the nearest-neighbor calcifications of two different clus- 
ters are within a cm in either the horizontal or vertical 
direction, the total cluster is counted as one; this is some- 
times defined as a diffuse cluster situation, and it admits the 
possibility of chaining clusters together. 

C. Tabulated Results 

The five sets of detection results are shown in Table 3. 
Each trial corresponds to a different threshold x or P/t). The 
thresholds corresponding to the five trials are arranged so 
that t 1 <t 2 < . . . <t 5 , and the corresponding values of PyxlO 4 
are in the last column. The goal is to identify the t where the 
experimental value of the sensitivity begins to drop below 
100%. In this table, the following definitions are used: 

Number of normals correctly classified 

Specificity = and 

Total number of normals 

Number of clusters found 
Sensitivity = 

Total number of abnormals 


TABLE 4 



THEORETICAL LIMITS FOR EACH TRIAL 

Trial 

Max 

Min 

Actual Counts 

P f (t) x 10 4 


63.92 

31.46 

44 

6.00 

*2 

31.46 

15.73 

31 

3.00 

% 

26.20 

13.10 

18 

2.50 

t 4 

10.48 

5.24 

13 

1.00 

% 

0.52 

0.26 

3 

0.05 


D. Observations 

The specificity rates (Table 3) are very encouraging, since 
a feasible operating Py can be found. In this case, it is trial 
3 and only one very subtle detection case is missed in trial 
4. This indicates we can hope to operate at 100% sensitivity 
while identifying 46% of the normal images. The theoretical 
isolated FP rates (Table 4) are in general agreement with the 
counted data. As the FP rate is reduced the agreement 
diverges somewhat because the integral required to find the 
FP rate is only an approximation. The final detection results 
are good order-of-magnitude estimates. These evaluation 
results indicate that the detection method behaves as 
predicted, and gives credence to the statistical modeling. If 
the model was merely a crude approximation, it is quite 
likely the detection results would not be in such close 
agreement. These detection results can be used to make 
inferences concerning processing full size images. For the 
isolated and cluster FP rate the conversion is 

Old rate X Full image breast area 

New rate = : 

Cut image area 


In going from trial 1 to trial 2, there is no measurable change 
in the evaluation. This means t was not changed enough. 
The parameters associated with trial 3 are the best, since it 
is possible to keep the sensitivity at 100%, and still identify 
46% of the normals. 

An estimate for the theoretical maximum number of 
isolated (individual) calcifications per image can be found 
by the formula 

18x256x256 + 12x256x128 
Max = 2 xP f x — 


and the minimum is given by: Min=Max/2. This formula 
comes from considering that there are two reduced images 
for each raw image and there are two possible sizes of 
reduced images (256x256) or (256x128). There are 18 large 
images and 12 small images. The Max is two times the Min 
because both d 3 and d 4 can contribute to detection, and it is 
possible to have no overlapping error in each image. These 
results are summarized in Table 4. 


TABLE 3 



EVALUATION OF EACH OF FIVE TRIALS 



Specificity 

Sensitivity 



Trial 

% 

% 

FP clusters/image 

Pf(T) X 10 4 

Ti 

15 

100 

1.20 

6.00 

**2 

15 

100 

1.36 

3.00 

^3 

46 

100 

0.93 

2.50 

t 4 

46 

94 

0.67 

1.00 

% 

92 

89 

0.13 

0.05 


Normally a comparison of detection results obtained by 
different techniques is useful. However, this is not appro- 
priate here for three reasons: (1) Using large image sections 
without artifacts is equivalent to processing ideal images. (2) 
The detection process was evaluated on a limited database 
consisting of 30 images. Excluding one or two images from 
the study may result in significant changes in detection 
results. (3) Very little other work has been done at 35 -/um 
resolution, and no work with common criteria that would 
allow a meaningful comparison. 

VIII. An Example Sequence of Steps for 
Application to an Image 

The following exemplary method can be repeated a 
plurality of times to find operating parameters using data 
that have known ground-truth information obtained from 
using the appropriate protocol. 

1. Locate the tissue region to be studied (e.g., a breast 
region) and eliminate the off-tissue region noise field and 
anomalous regions. 

2. Apply a wavelet expansion and separate at least one 
important component expansion image (for the example of 
breast microcalcifications, two expansion images subimages 
are separated out). 

3. Determine the functional probability form for each 
expansion image. 

4. Using maximum likelihood arguments, determine the 
optimal statistical normal image region test for each expan- 
sion image. 

5. Apply the test to a training data set using an initial 
estimate for the thresholds, where in the most general terms, 
the operating threshold is different in each expansion image. 
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a. Test each region in the first expansion image (e.g., each 
8x8 region in the d 3 image). If the calculated test 
statistic is less than the threshold, set the region to zero 
because it is assumed normal; else set it to one. 

b. Apply the same test to each successive expansion image 
selected (e.g., each 16x16 region in the d 4 ) image. 

6. Form the union of steps 5a and 5b as the total output 
binary mask. 

7. Compare the detection results for all images with the 
ground-truth data and repeat the process until optimal oper- 
ating parameters are found. 

8. Once the parameters are found, the approach can be 
applied to the test set following steps 1-6. The process can 
terminate here or continue. 

9. Following step 6, a simple counting procedure can be 
implemented if desired to eliminate isolated suspicious 
areas, since fewer than three or five suspicious areas within 
1 cm 2 proximity are assumed to be not important. 

10. Project the binary detection mask into the raw image 
or any combination of detail images for finer classification 
algorithms (see FIGS. 17-23). 

Alternatively, the method can be used as a calcification 
detection algorithm by adjusting the detection operating 
parameters (step 5) with different optimization goals. 

IX. Conclusions 

The detection technique presented here is dictated by the 
data and follows from the Neyman-Pe arson lemma. If the 
pdf model is a reasonable approximation, this is the most 
powerful test for finding regions that deviate significantly 
from normal regions, that is, regions that may contain 
calcifications or other artifacts of similar scale. In this sense, 
the test is optimal. In essence, the technique merges two 
powerful analysis techniques: classical signal detection 
theory and multiresolution decomposition. 

In general, it is illustrated herein that mammograms are 
not as irregular as commonly believed. The evidence clearly 
indicates that these mammograms share common statistical 
qualities. In the vast majority of images studied, a primary 
and summary statistics appear to be dependable estimators 
for the detection scheme. 

This method of differentiating normal regions from poten- 
tially abnormal regions containing calcifications should be 
considered as the initial stage of processing. It is important 
that this separation technique be optimal; if it fails, any 
ensuing processing will naturally fail. As illustrated by the 
detection analysis, in order to identify a significant fraction 
of the normal images while maintaining a low FN rate 
(essentially zero), the isolated FP error needs to be reduced. 
One way to address this error is with two additional stages 
of processing: (al) examine each small suspicious region in 
more detail for further discrimination, and (a2) consider 
regions on the order of a square centimeter and eliminate the 
entire region if the total suspicious area within the region 
does not meet the cluster criterion of three classifications 
within a square centimeter. In order to overcome intrinsic FP 
error (see the summary statistic), the decision criteria for the 
first additional stage should not be strongly dependent on the 
previous decision mechanism. The detection mask approach 
renders itself to further analysis of type (al), because more 
information can be included in the suspicious regions by 
projecting the mask into a more resolved image, as illus- 
trated in FIG. 7. The second additional stage, (a2), can be 
implemented by a simple counting procedure. These addi- 
tional stages are being developed. 

The detection process was illustrated with the symmlet 
basis. Other wavelet bases can be used for comparison 
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purposes to optimize the choice of bases. Thresholds can be 
set the same and the experiment repeated. The 30 images are 
a fair representation of a clinical mammography database. 
Thus using the present invention, one can quantify the “best 
5 basis” for usage in mammography. 

The evaluation results provide a strong impetus for further 
pursuit and analysis of the multiresolution statistical tech- 
nique. An automated method should be developed to seg- 
regate the isolated calcifications that do not belong to a 
io cluster. 

It may be appreciated by one skilled in the art that 
additional embodiments may be contemplated, including in 
the analysis of other types of medical (and other) images. 

In the foregoing description, certain terms have been used 
15 for brevity, clarity, and understanding, but no unnecessary 
limitations are to be implied therefrom beyond the require- 
ments of the prior art, because such words are used for 
description purposes herein and are intended to be broadly 
construed. Moreover, the embodiments of the apparatus 
20 illustrated and described herein are by way of example, and 
the scope of the invention is not limited to the exact details 
of construction. 

Having now described the invention, the construction, the 
operation and use of preferred embodiment thereof, and the 
25 advantageous new and useful results obtained thereby, the 
new and useful constructions, and reasonable mechanical 
equivalents thereof obvious to those skilled in the art, are set 
forth in the appended claims. 
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What is claimed is: 

1. A method for analyzing a medical image to determine 
whether the image is classifiable as normal, the method 
comprising the steps of: 

applying a wavelet expansion to a digital representation of 
a raw image, the raw image comprising an array of 
sectors, each sector having an intensity level, to obtain 
a plurality of expansion images of varying resolution; 

selecting at least one expansion image having a resolution 
commensurate with a desired predetermined detection 
resolution range; 

dividing each expansion image into a plurality of regions, 
each region comprising at least one sector; and 

creating an output image comprising a combination of all 
regions for each selected expansion image, each region 
having a first value when the region intensity level is 
above a predetermined threshold level and a second 
value when the region intensity level is below the 
threshold level, for localizing a potential abnormality 
within the image; 

wherein an absence of a predetermined number of regions 
having a first value intensity level is indicative of the 
image being classifiable as normal. 

2. The method recited in claim 1, further comprising the 
steps, prior to the output image creating step, of: 

applying a wavelet expansion to a digital representation of 
a raw control image comprising an array of sectors, 
each sector having an intensity level, to obtain a 
plurality of expansion images of varying resolution; 

selecting at least one expansion image having a resolution 
commensurate with a desired predetermined detection 
resolution range; 


determining a functional form of a probability distribution 
function for each selected expansion image; 

determining an optimal statistical normal image region 
test for each selected expansion image; 

5 establishing the threshold level for the probability distri- 
bution function from the optimal statistical normal 
image region test for each selected expansion image. 

3. The method recited in claim 2, wherein the expansion- 
image selecting step comprises empirically selecting two 

10 adjacent expansion images. 

4. The method recited in claim 3, wherein the test- 
determining step comprises forming a search window hav- 
ing a predetermined resolution size for each selected expan- 
sion image, and wherein the output image creating step 

15 comprises combining the selected two adjacent expansion 
images and applying a binary mask having the predeter- 
mined resolution size to the combined expansion space 
images. 

5. The method recited in claim 2, wherein the optimal 

20 statistical normal image region test determining step com- 
prises using a maximum likelihood technique. 

6. The method recited in claim 2, wherein the threshold 
level establishing step comprises applying an iterative pro- 
cedure to determine a set of operating parameters for each 

25 selected expansion image. 

7. The method recited in claim 2, wherein the probability 
density function comprises one of a family of distributions 
having the form: 



8. The method recited in claim 2, wherein the test deter- 

35 mining step comprises assuming a parametric form of a test 

statistic distribution. 

9. The method recited in claim 1, wherein the test deter- 
mining step comprises assuming a nonparametric form of a 
test statistic distribution and using a kernel density estima- 

40 tion method. 

10. The method recited in claim 1, wherein the raw image 
comprises a film image, and further comprising the step, 
prior to the wavelet-expansion application step, of digitizing 
the film image. 

45 11. The method recited in claim 1, wherein the wavelet- 

expansion application step comprises performing a sepa- 
rable kernel two-dimensional pyramid downsampling/ 
upsampling scheme. 

12. The method recited in claim 11, wherein each expan- 

50 sion image comprises a difference in information between a 

first image having a first resolution and a second image 
having a resolution one-half that of the first image. 

13. The method recited in claim 1, wherein the expansion- 
image selecting step comprises empirically selecting a num- 

55 ber of expansion images commensurate with an initial raw 
image resolution. 

14. The method recited in claim 1, wherein the expansion 
image selecting step comprises selecting a plurality of 
expansion images, each expansion image having a resolu- 

60 tion commensurate with a size range of an abnormality 
desired to be detected. 

15. The method recited in claim 1, further comprising the 
step of defining a tissue boundary for analysis. 

16. The method recited in claim 15, wherein the medical 

65 image comprises a digitized mammogram, and wherein the 

tissue boundary defining step comprises excising a breast 
boundary on the digitized mammogram. 
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17. The method recited in claim 16, wherein the breast 
boundary excising step comprises the steps of: 

defining an off-tissue region comprising a generally ran- 
dom noise field and a plurality of anomalous regions; 

defining a tissue region comprising signal information; 

separating out the random noise field; 

determining a remaining contiguous region containing 
information; and 

setting the image region outside the remaining contiguous 
region to zero. 

18. The method recited in claim 1, further comprising the 
steps of: 

defining a cluster area having a size greater than a size of 
the region; 

setting a number of regions within each cluster area in the 
image having the first value to a first variable; 

comparing the first variable with a predetermined second 
variable; 

if the first variable is greater than or equal to the second 
variable, flagging the cluster area as potentially suspi- 
cious for the presence of an abnormality; and 

if the first variable is less than the second variable, 
classifying the cluster area to be normal. 

19. The method recited in claim 1, wherein the medical 
image comprises a digitized mammogram, and the expan- 
sion image selecting step comprises selecting two adjacent 
expansion images having resolution less than 0.5 mm. 

20. The method recited in claim 19, wherein the two 
selected adjacent expansion images further have resolution 
greater than 0.1 mm. 

21. An apparatus for analyzing a medical image to deter- 
mine whether an abnormality is present comprising: 

means for applying a wavelet expansion to a digital 
representation of a raw image, the raw image compris- 
ing an array of sectors, each sector having an intensity 
level, to obtain a plurality of expansion images of 
varying resolution; 

means for selecting at least one expansion image having 
a resolution commensurate with a desired predeter- 
mined detection resolution range; 
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means for dividing each expansion image into a plurality 
of regions, each region comprising at least one sector; 
and 

means for creating an output image comprising a combi- 
nation of all regions for each selected expansion image, 
each region having a first value when the region 
intensity level is above a predetermined threshold level 
and a second value when the region intensity level is 
io below the threshold level, for localizing a potential 
abnormality within the image. 

22. The apparatus recited in claim 21, further comprising: 
means for applying a wavelet expansion to a digital 

representation of a raw control image, the raw control 
15 image comprising an array of sectors, each sector 
having an intensity level, to obtain a plurality of 
expansion images of varying resolution; 
means for selecting at least one expansion image having 
20 a resolution commensurate with a desired predeter- 
mined detection resolution range; 
means for determining a functional form of a probability 
distribution function for each selected expansion 
image; 

25 

means for determining an optimal statistical normal image 
region test for each selected expansion image; and 
means for establishing the threshold level for the prob- 
ability distribution function from the optimal statistical 
30 normal image region test for each selected expansion 
image. 

23. The apparatus recited in claim 21, wherein the raw 
image comprises a film image, and further comprising 
means for digitizing the raw image. 

35 24. The apparatus recited in claim 21, wherein the wavelet 

expansion applying means comprises software means. 

25. The apparatus recited in claim 21, wherein the output 
image creating means comprises software means resident 
within a processor and a screen in electronic communication 
40 with the processor for visualizing the output image. 



